% Constructing synthetic series for nonmarket work and leisure, based on BLS
% actual hours worked per week at the respondent's main job
% for our local projection regressions with uncertainty shocks
% Daniela Hauser, February 2022
% Updated by James Cabral, March 2023 to use aggregate hours: nonfarm
% payrolls

clear all
close all
clc

%addpath("MatlabFunctions")
addpath("../data/synthetic_series")
addpath("../data/ATUS_AHTUS_state_datasets")

%% 1) Required inputs and load data

% Set the specific quarter of ATUS nonmarket hours/worked per week (imported as `Data_ATUS_Raw' below) 
% to be used to recover levels of nonmarket hours worked

ATUSobs = 68; % last available ATUS observation 2019Q4
BLSobs = 236; %observation in the BLS data corresponding to 2019Q4
startyear = 1960

% Choose the following two options for final table summarizing business-cycle moments:
% 1) choose detrending method 
detrendM = 1; % linear trend

%Generate date arrays which will be used for the YoY and QoQ series
%These will adapt the code for the years in which we want the synthetic
%series

DateStringQuart_1 = {};

for yr = startyear+1:2021
    if yr == 2021
        DateStringQuart_1{end+1,1} = [num2str(yr), 'Q', num2str(1)];
    else
        for Q = 1:4
                DateStringQuart_1{end+1,1} = [num2str(yr), 'Q', num2str(Q)];
        end
    end
end

DateStringQuart_2 = {};
 
for yr = startyear:2021    
    if yr == 2021
        DateStringQuart_2{end+1,1} = [num2str(yr), 'Q', num2str(1)]
    else
        if yr == startyear 
            for Q = 2:4
                DateStringQuart_2{end+1,1} = [num2str(yr), 'Q', num2str(Q)];
            end
        else
            for Q = 1:4
                DateStringQuart_2{end+1,1} = [num2str(yr), 'Q', num2str(Q)];
            end
        end
    end
end

formatInQuart   = 'yyyyQQ';
dates_Quart_1   = datenum(DateStringQuart_1,formatInQuart);
dates_Quart_2   = datenum(DateStringQuart_2,formatInQuart);           

% Load raw data series of per capita nonfarm market hours worked
BLS_Raw       = xlsread('NONFARM_HOURS_Q.xls','Sheet1');
BLS_Raw       = BLS_Raw(startyear-1959:245,3);
BLS_Raw_Level = BLS_Raw;
BLS_Raw       = log(BLS_Raw);

% Estimated substitutability of nonmarket work and leisure
% based on regression specification as in Aguiar et al. 2013
% using bi-annual averages, for sample 2003-2019
% see Table 1 in our working document
BETA_Nonmarket  = -0.2530; %nonmarket work
BETA_Leisure    = -0.6001; %leisure

% Load aggregate ATUS series for nonmarket hours worked and leisure, per week (Raw Data)
Data_ATUS_Raw = xlsread('state_dataset_ALLUS.xls');
Data_ATUS_Raw = [Data_ATUS_Raw(:,30) Data_ATUS_Raw(:,38)];

%% 3) Compute differences of market hours worked for different detrending methods
% Compute quarterly year-over-year differences of market hours worked
BLS_Lead4_Level   =   BLS_Raw_Level(5:end,:);
D_BLS_YoverY_Level=   BLS_Lead4_Level- BLS_Raw_Level(1:end-4,:);


%% 4) Compute the aggregate nonmarket hours worked (cyclical variation)
%     by multiplying the estimated substituability (BETA_Nonmarket) with the aggregate changes in market hours worked

% Year-over-year changes at quarterly differences
D_Nonmarket_YoverY_Level    = BETA_Nonmarket*D_BLS_YoverY_Level(:,:);
D_Leisure_YoverY_Level      = BETA_Leisure*D_BLS_YoverY_Level(:,:);

Annual_Betas =readtable('ATUS_AHTUS_betas.csv');

%Create a row of betas for nonmarket work and leisure to be multiplied with the  
BETA_Nonmarket_Annual = cell(28, 2);
BETA_Leisure_Annual = cell(28, 2);

%set up the beta series so that each year is matched to the correct beta
for i = startyear:2021
    if (i >= startyear) & (i <1998)
        BETA_Nonmarket_Annual{i-(startyear-1), 1} = table2array(Annual_Betas((Annual_Betas.timeuse == "Nonmarket Work" & Annual_Betas.year == 1995), 3))
        BETA_Leisure_Annual{i-(startyear-1), 1} = table2array(Annual_Betas((Annual_Betas.timeuse == "Leisure" & Annual_Betas.year == 1995), 3))
        BETA_Nonmarket_Annual{i-(startyear-1), 2} = i
        BETA_Leisure_Annual{i-(startyear-1), 2} = i
        
    elseif (i >= 1998) & (i <2001)
        BETA_Nonmarket_Annual{i-(startyear-1), 1} = table2array(Annual_Betas((Annual_Betas.timeuse == "Nonmarket Work" & Annual_Betas.year == 1998), 3))
        BETA_Leisure_Annual{i-(startyear-1), 1} = table2array(Annual_Betas((Annual_Betas.timeuse == "Leisure" & Annual_Betas.year == 1998), 3))
        BETA_Nonmarket_Annual{i-(startyear-1), 2} = i
        BETA_Leisure_Annual{i-(startyear-1), 2} = i 
        
    elseif (i >= 2001) & (i < 2004)
        BETA_Nonmarket_Annual{i-(startyear-1), 1} = table2array(Annual_Betas((Annual_Betas.timeuse == "Nonmarket Work" & Annual_Betas.year == 2004), 3))
        BETA_Leisure_Annual{i-(startyear-1), 1} = table2array(Annual_Betas((Annual_Betas.timeuse == "Leisure" & Annual_Betas.year == 2004), 3))        
        BETA_Nonmarket_Annual{i-(startyear-1), 2} = i
        BETA_Leisure_Annual{i-(startyear-1), 2} = i    
        
    elseif (i >= 2004) & (i <=2019)
        BETA_Nonmarket_Annual{i-(startyear-1), 1} = table2array(Annual_Betas((Annual_Betas.timeuse == "Nonmarket Work" & Annual_Betas.year == i), 3))
        BETA_Leisure_Annual{i-(startyear-1), 1} = table2array(Annual_Betas((Annual_Betas.timeuse == "Leisure" & Annual_Betas.year == i), 3)) 
        BETA_Nonmarket_Annual{i-(startyear-1), 2} = i
        BETA_Leisure_Annual{i-(startyear-1), 2} = i
    elseif (i > 2019)
        BETA_Nonmarket_Annual{i-(startyear-1), 1} = table2array(Annual_Betas((Annual_Betas.timeuse == "Nonmarket Work" & Annual_Betas.year == 2019), 3))
        BETA_Leisure_Annual{i-(startyear-1), 1} = table2array(Annual_Betas((Annual_Betas.timeuse == "Leisure" & Annual_Betas.year == 2019), 3)) 
        BETA_Nonmarket_Annual{i-(startyear-1), 2} = i
        BETA_Leisure_Annual{i-(startyear-1), 2} = i
    end
   
end

%set of betas for annual datasets
BETA_Nonmarket_A_Annual = BETA_Nonmarket_Annual(2:end,:)
BETA_Leisure_A_Annual = BETA_Leisure_Annual(2:end,:)

%set up consistent variable names
BETA_Nonmarket_A_Annual_Yrs = cell2table(BETA_Nonmarket_A_Annual, 'VariableNames', ["beta", "Year"]);
BETA_Leisure_A_Annual_Yrs = cell2table(BETA_Leisure_A_Annual, 'VariableNames', ["beta", "Year"]);

D_BLS_YoverY_Level_Yrs = table(D_BLS_YoverY_Level, str2double(extractBetween(DateStringQuart_1, 1, 4)), 'VariableNames', ["MktHrs", "Year"]);

%now inner join by year and multiply these betas by the YoY differences
%to reverse transformations from Stata output, we must multiply by -1 and
%divide by 100
T =  innerjoin(D_BLS_YoverY_Level_Yrs, BETA_Nonmarket_A_Annual_Yrs);
D_Nonmarket_YoverY_Level_Annual = table2array(T(:,1)).*(table2array(T(:,3))/100*(-1));

TT =  innerjoin(D_BLS_YoverY_Level_Yrs, BETA_Leisure_A_Annual_Yrs);
D_Leisure_YoverY_Level_Annual = table2array(TT(:,1)).*(table2array(TT(:,3))/100*(-1));

clear T TT; 

%% 5) Converting differences to levels - Year-over-year differences

% Set initial values used to convert year-over-year changes to levels
% Here we need 4 initial values. 
% We therefore use the same observations as for year-over-year changes 
% (the one chosen in ATUSobs) for all four observations
Nonmarket_YoverY_Level(BLSobs-3,:)              = Data_ATUS_Raw(ATUSobs,2); 
Nonmarket_YoverY_Level(BLSobs-2,:)              = Data_ATUS_Raw(ATUSobs,2); 
Nonmarket_YoverY_Level(BLSobs-1,:)              = Data_ATUS_Raw(ATUSobs,2); 
Nonmarket_YoverY_Level(BLSobs  ,:)              = Data_ATUS_Raw(ATUSobs,2); 

Leisure_YoverY_Level(BLSobs-3,:)                = Data_ATUS_Raw(ATUSobs,1); 
Leisure_YoverY_Level(BLSobs-2,:)                = Data_ATUS_Raw(ATUSobs,1); 
Leisure_YoverY_Level(BLSobs-1,:)                = Data_ATUS_Raw(ATUSobs,1); 
Leisure_YoverY_Level(BLSobs  ,:)                = Data_ATUS_Raw(ATUSobs,1); 

%set same initial values for series using annual betas
%once again, annual betas are only used for nonmarket hours and leisure
Nonmarket_YoverY_Level_Annual(BLSobs-3,:)              = Data_ATUS_Raw(ATUSobs,2); 
Nonmarket_YoverY_Level_Annual(BLSobs-2,:)              = Data_ATUS_Raw(ATUSobs,2); 
Nonmarket_YoverY_Level_Annual(BLSobs-1,:)              = Data_ATUS_Raw(ATUSobs,2); 
Nonmarket_YoverY_Level_Annual(BLSobs  ,:)              = Data_ATUS_Raw(ATUSobs,2); 

Leisure_YoverY_Level_Annual(BLSobs-3,:)                = Data_ATUS_Raw(ATUSobs,1); 
Leisure_YoverY_Level_Annual(BLSobs-2,:)                = Data_ATUS_Raw(ATUSobs,1); 
Leisure_YoverY_Level_Annual(BLSobs-1,:)                = Data_ATUS_Raw(ATUSobs,1);
Leisure_YoverY_Level_Annual(BLSobs  ,:)                = Data_ATUS_Raw(ATUSobs,1); 


% Moving backwards - year-over-year differences
for jjj = 4:BLSobs-1
    Nonmarket_YoverY_Level(BLSobs-jjj,:)        = Nonmarket_YoverY_Level(BLSobs+4-jjj,:)  - D_Nonmarket_YoverY_Level(BLSobs+4-jjj,:);
    Leisure_YoverY_Level(BLSobs-jjj,:)          = Leisure_YoverY_Level(BLSobs+4-jjj,:)    - D_Leisure_YoverY_Level(BLSobs+4-jjj,:);
    Nonmarket_YoverY_Level_Annual(BLSobs-jjj,:)        = Nonmarket_YoverY_Level_Annual(BLSobs+4-jjj,:)  - D_Nonmarket_YoverY_Level_Annual(BLSobs+4-jjj,:);
    Leisure_YoverY_Level_Annual(BLSobs-jjj,:)          = Leisure_YoverY_Level_Annual(BLSobs+4-jjj,:)    - D_Leisure_YoverY_Level_Annual(BLSobs+4-jjj,:);    
end

% Moving forwards
for iii = BLSobs+1:size(D_Nonmarket_YoverY_Level,1)
    Nonmarket_YoverY_Level(iii,:)        = D_Nonmarket_YoverY_Level(iii,:)        + Nonmarket_YoverY_Level(iii-3,:);
    Leisure_YoverY_Level(iii,:)          = D_Leisure_YoverY_Level(iii,:)          + Leisure_YoverY_Level(iii-3,:);

end

% Moving forwards (using annual betas)
for iii = BLSobs+1:size(D_Nonmarket_YoverY_Level_Annual,1)
    Nonmarket_YoverY_Level_Annual(iii,:)        = D_Nonmarket_YoverY_Level_Annual(iii,:)        + Nonmarket_YoverY_Level_Annual(iii-3,:);
    Leisure_YoverY_Level_Annual(iii,:)          = D_Leisure_YoverY_Level_Annual(iii,:)          + Leisure_YoverY_Level_Annual(iii-3,:);
end

writetable(table(Nonmarket_YoverY_Level(:,1),'VariableNames', ["nonmkt"]), '../data/synthetic_series/nonmkt_YoY_q_nonfarm.csv')
writetable(table(Leisure_YoverY_Level(:,1),'VariableNames', ["leisure"]), '../data/synthetic_series/leisure_YoY_q_nonfarm.csv')
writetable(table(Nonmarket_YoverY_Level_Annual(:,1),'VariableNames', ["nonmkt_an"]), '../data/synthetic_series/nonmkt_YoY_q_annual_beta_nonfarm.csv')
writetable(table(Leisure_YoverY_Level_Annual(:,1),'VariableNames', ["leisure_an"]), '../data/synthetic_series/leisure_YoY_q_annual_beta_nonfarm.csv')

disp('DONE')